#include "slepbase.h"
#include "epph.h"


void glLogisticR(slep_matrix* A, slep_vector* y, double z, OPTS* opts,
				 slep_vector* x, double *c, slep_vector* funVal, slep_vector* valueL){

	double sWeight[2], m1=0., m2=0., k1,k2, lambda, lambda_max;
	slep_vector *weight, *b, *x1, *Ax, *ATb, *b1, temp;
	int *ind, ind_size;
	int i,k;
	double q;
	slep_vector* gWeights;

	assert( y->dim == A->m);

	initNormalization(A,opts);

	weight = slepAllocDenseVector(y->dim);
	b = slepAllocDenseVector(y->dim);
	x1 = slepAllocDenseVector(x->dim);
	b1 = slepAllocDenseVector(y->dim);
	ATb = slepAllocDenseVector(x->dim);
	Ax = slepAllocDenseVector(y->dim);

	if( ! hasGot(opts,SlepOPT_indices)){
		assert( (opts->REGISTER & SlepOPT_indices)== SlepOPT_indices);
	}else {
		ind = opts->ind;
		ind_size = opts->ind_size;
		k = ind_size - 1;
		assert( ind[k] == A->n);
	}

	if( ! hasGot(opts,SlepOPT_q))
		q=2.0, opts->q = 2.0, opts->REGISTER = ((opts->REGISTER)|SlepOPT_q);
	else{
		q = opts->q;
		assert( q >= 1.0);
	}

	if(hasGot(opts,SlepOPT_sWeight)){
		sWeight[0] = opts->sWeight[0];
		sWeight[1] = opts->sWeight[1];
		for(i=0; i < y->dim; i++)
			m1+= ((y->data.vctr[i])>0?1.:0.);
		m2 = i-m1;
		m1*=sWeight[0];
		m2*=sWeight[1];

		k1=sWeight[0]/(m1+m2);
		k2=sWeight[1]/(m1+m2);
		weight = slepAllocDenseVector(y->dim);
		m1=m2=0.0;
		for(i=0; i < weight->dim; i++){
			(weight->data.vctr[i]) = ((y->data.vctr[i])>0?k1:k2);
			m1 += ((y->data.vctr[i])>0?(weight->data.vctr[i]):0.0);
			m2 += ((y->data.vctr[i])>0?0.0:(weight->data.vctr[i]));
		}
	}else{
		weight = slepAllocDenseVector(y->dim);
		m1= 1.0 / ((double)(y->dim));
		for(i=0; i < y->dim; i++){
			m2+= ((y->data.vctr[i])>0?0.:1.);
			(weight->data.vctr[i]) = m1;
		}
		m1 = i-m2;
		m1 *= (1.0 / ((double)(y->dim)));
		m2 *= (1.0 / ((double)(y->dim)));
	}

	if( hasGot(opts, SlepOPT_gWeights)){
		gWeights = opts->gWeights;
		assert( gWeights->dim == k);
		slep_dscal(-1.0,gWeights);
		assert( slep_dmax(gWeights) <= 0.0 );
		slep_dscal(-1.0,gWeights);
	}else {
		gWeights = slepAllocDenseVector(k);
		DMEMSET(gWeights->data.vctr,1.0,k,i);
	}

	if( opts->rFlag == 0)
		lambda = z;
	else{
		slep_vector* norm_ATb;
		double q_bar;
		for(i=0; i < y->dim; i++)
			b->data.vctr[i] = ((y->data.vctr[i])>0 ? m2 : -m1);
		slep_xdiagy(weight,b);

		NORMALIZEDMTV(A,b,ATb,opts, b1);	


		if( q==1.0)
			q_bar = 1e30;
		else if( q >= 1e6)
			q_bar = 1.0;
		else
			q_bar = q/(q-1.0);
		norm_ATb = slepAllocDenseVector(k);
		slep_zerov(norm_ATb);
		for(i=0; i < k; i++){
			temp = (*ATb);
			temp.data.vctr = ATb->data.vctr +(ind[i]);
			temp.dim = ind[i+1] - ind[i];
			norm_ATb->data.vctr[i] = slep_dpnrm(&temp,q_bar);
		}

		slep_xdiagy(gWeights,norm_ATb);
		lambda_max = slep_dmax(norm_ATb);
		lambda = z * lambda_max;
		slepFreeVector(norm_ATb);
	}

	if( opts->init == 2){
		slep_zerov(x);
		(*c) = log( m1/m2);
	}else{
		if(hasGot(opts,SlepOPT_x0))
			slep_dvcopy(opts->x0, x);
		else
			slep_zerov(x);
		if(hasGot(opts, SlepOPT_c0))
			(*c) = ((opts->c0)->data.vctr[0]);
		else
			(*c) = log( m1/m2);
	}

	NORMALIZEDMV(A,x,Ax,opts,x1,i);

	/* The main program */ 

	if( opts->mFlag == 0 && opts->lFlag==0){
		int bFlag=0;
		slep_vector* rho;
		double L = 1.0/(A->m);
		double cp, ccp, alphap, alpha, beta, sc;
		slep_vector *weighty, *Axp, *xp, *xxp;
		slep_vector *s, *aa, *As, *bb, *prob, *g, *v, *norm_x_k;
		int iterStep;
		double fun_s, gc, l_sum, r_sum, fun_x;
		
		weighty = slepCopyToNewVector(weight);
		slep_xdiagy(y,weighty);

		xp = slepCopyToNewVector(x);
		Axp = slepCopyToNewVector(Ax);
		xxp = slepAllocDenseVector(x->dim);
		slep_zerov(xxp);
		cp = (*c);
		ccp = 0;

		/* The Armijo Goldstein line search scheme + accelearted gradient descent */

		alphap = 0.; alpha = 1.;

		s = slepAllocDenseVector(x->dim);
		As = slepAllocDenseVector(Ax->dim);
		aa = slepAllocDenseVector(As->dim);
		bb = slepAllocDenseVector(As->dim);
		prob = slepAllocDenseVector(As->dim);
		g = slepAllocDenseVector(ATb->dim);
		v = slepAllocDenseVector(x->dim);
		rho = slepAllocDenseVector(gWeights->dim);
		norm_x_k = slepAllocDenseVector(k);

		for(iterStep = 0; iterStep< opts->maxIter; iterStep ++){
			beta = (alphap - 1.)/ alpha;
			slep_dvcopy(x,s);
			slep_daxpy(beta,xxp,s);
			sc = (*c) + beta * ccp;

			slep_dvcopy(Ax,As);
			slep_daxpby(-beta,Axp,1.0+beta, As);
			
			slep_dvcopy(As, aa);
			slep_xdiagy(y,aa);
			slep_daxpby(-sc,y,-1.0,aa);

			gc = 0.;
			fun_s = 0.;
			for(i=0; i< aa->dim; i++){
				bb->data.vctr[i] = max(aa->data.vctr[i], 0.);
				fun_s +=  weight->data.vctr[i] * ( log( 
														exp(-(bb->data.vctr[i]))+
														exp((aa->data.vctr[i])-(bb->data.vctr[i])))
													+ (bb->data.vctr[i])) ;
				prob->data.vctr[i] = 1. / (1. + exp(aa->data.vctr[i]) );
				b->data.vctr[i] = -(weighty->data.vctr[i]) *(1 - prob->data.vctr[i]);

				gc += b->data.vctr[i];
			}

			NORMALIZEDMTV(A,b,g,opts,b1);

			slep_dvcopy(x,xp);
			slep_dvcopy(Ax,Axp);
			cp = (*c);

			for(;;){
				slep_dvcopy(s,v);
				slep_daxpy(-1.0/L, g, v);

				(*c) = sc - gc/L;

				slep_dvcopy(gWeights,rho);
				slep_dscal(lambda/L, rho);
				if( q< 1e6)
					eppVectorI(x->data.vctr,v->data.vctr, ind,k, x->dim,rho->data.vctr, q);
				else
					eppVectorI(x->data.vctr,v->data.vctr, ind,k, x->dim,rho->data.vctr, 1.e6);

				for(i=0; i < x->dim; i++)
					v->data.vctr[i] = x->data.vctr[i] - s->data.vctr[i];

				NORMALIZEDMV(A,x,Ax,opts,x1,i);

				slep_dvcopy(Ax, aa);
				slep_xdiagy(y,aa);
				slep_daxpby(-(*c),y,-1.0,aa);

				fun_x = 0.;
				for(i=0; i< aa->dim; i++){
					bb->data.vctr[i] = max(aa->data.vctr[i], 0.);
					fun_x +=  weight->data.vctr[i] * ( log( 
															exp(-(bb->data.vctr[i]))+
															exp((aa->data.vctr[i])-(bb->data.vctr[i])))
														+ (bb->data.vctr[i])) ;
				}

				r_sum = (slep_ddot(v,v) + ((*c)-sc)*((*c)-sc))/2.0;
				l_sum = fun_x - fun_s - slep_ddot(v,g) - ((*c)-sc)* gc ;

				if( r_sum <= 1e-20 ){
					bFlag=1;
					break;
				}

				if( l_sum <= r_sum * L)
					break;
				else
					L = max(2.0*L, l_sum/r_sum);
			}

			alphap = alpha;
			alpha = (1.0 +sqrt(4.0*alpha * alpha + 1))/2.0;

			valueL->data.vctr[iterStep] = L;
			slep_dvcopy(x,xxp);
			slep_daxpy(-1.0, xp, xxp);
			ccp = (*c) - cp;

			slep_zerov(norm_x_k);

			for(i=0; i < k; i++){
				temp = (*x);
				temp.data.vctr = x->data.vctr +(ind[i]);
				temp.dim = ind[i+1] - ind[i];
				norm_x_k->data.vctr[i] = slep_dpnrm(&temp,q);
			}
			funVal->data.vctr[iterStep] = fun_x + lambda * slep_ddot(norm_x_k,gWeights);

			if(bFlag)
				break;

			if( opts->tFlag == 0){
				if( iterStep >= 2){
					if( fabs(funVal->data.vctr[iterStep] - funVal->data.vctr[iterStep-1]) < opts->tol)
						break;
				}
			}else if( opts->tFlag == 1){
				if(iterStep>=2){
					if( fabs(funVal->data.vctr[iterStep] - funVal->data.vctr[iterStep-1]) < (opts->tol)*(funVal->data.vctr[iterStep-1]) )
						break;
				}
			}else if(opts->tFlag == 2){
				if( funVal->data.vctr[iterStep] <= opts->tol)
					break;
			}else if(opts->tFlag == 3){
				if(slep_dnrm2(xxp) <= opts->tol)
					break;
			}else if(opts->tFlag == 4){
				if( slep_dnrm2(xxp) <= (opts->tol) * max(slep_dnrm2(xp),1))
					break;
			}else if(opts->tFlag == 5){
				if( iterStep >= opts->maxIter)
					break;
			}
		}
		slepFreeVector(weighty);
		slepFreeVector(Axp);
		slepFreeVector(xp);
		slepFreeVector(xxp);
		slepFreeVector(s);
		slepFreeVector(aa);
		slepFreeVector(As);
		slepFreeVector(bb);
		slepFreeVector(prob);
		slepFreeVector(g);
		slepFreeVector(v);
		slepFreeVector(norm_x_k);

		funVal->dim = iterStep;
		valueL->dim = iterStep;
	}

	/* Reformulated problem + Nemirovski's line search scheme */
	/*
	 * Original problem:           min  f(x,c) = - weight_i * log (p_i) + 1/2 rsL2 * ||x||_2^2 
	 *                                              + z * \|x\|_1
	 *
	 * Reformulated problem:   min  f(x,c) = - weight_i * log (p_i) + 1/2 rsL2 * ||x||_2^2
	 *                                                   + z * 1^T t
	 *                                   s.t., |x| <=t
	 **/

	if( opts->mFlag==1 && opts->lFlag==0  && (opts->q > 1.99 && opts->q < 2.01) ){
		int bFlag = 0, iterStep;
		double L = 1.0/(A->m);
		slep_vector *weighty, *xp, *Axp, *xxp, *t, *tp;
		double cp, ccp, alpha, alphap, beta;
		slep_vector *s, *As, *s_t, *aa, *bb, *g, *prob, *u, *v;
		double fun_s, fun_x, gc, sc, l_sum, r_sum;
		slep_vector *v_t;

		weighty = slepCopyToNewVector(weight);
		slep_xdiagy(y,weighty);
		xp = slepCopyToNewVector(x);
		Axp = slepCopyToNewVector(Ax);
		xxp = slepAllocDenseVector(x->dim);
		slep_zerov(xxp);
		cp = (*c);
		ccp = 0.;

		t = slepAllocDenseVector(k);
		slep_zerov(t);
		for(i=0; i < k; i++){
			temp = (*x) ; 
			temp.data.vctr = x->data.vctr +(ind[i]);
			temp.dim = ind[i+1] - ind[i];
			t->data.vctr[i] = slep_dnrm2(&temp);
		}
		tp = slepCopyToNewVector(t);

		alphap = 0. ;
		alpha = 1.;

		s = slepAllocDenseVector(x->dim);
		s_t = slepAllocDenseVector(t->dim);
		v = slepAllocDenseVector(s_t->dim);
		v_t = slepAllocDenseVector(t->dim);
		As = slepAllocDenseVector(Ax->dim);
		aa = slepAllocDenseVector(Ax->dim);
		bb = slepAllocDenseVector(Ax->dim);
		prob = slepAllocDenseVector(Ax->dim);
		g = slepAllocDenseVector(s->dim);
		u = slepAllocDenseVector(s->dim);

				
		for( iterStep = 0; iterStep < opts->maxIter; iterStep++){
			beta = (alphap - 1.) / alpha;

			slep_dvcopy(x,s);
			slep_daxpy(beta, xxp, s);
			sc = (*c) + beta * ccp;
			
			slep_dvcopy(t, s_t);
			slep_daxpby(-beta, tp,1.+ beta,s_t);

			slep_dvcopy(Ax,As);
			slep_daxpby(-beta, Axp,1.+beta, As);

			slep_dvcopy(As, aa);
			slep_xdiagy(y,aa);
			slep_daxpby(-sc,y,-1.0,aa);

			gc = 0.;
			fun_s = 0.;
			for(i=0; i< aa->dim; i++){
				bb->data.vctr[i] = max(aa->data.vctr[i], 0.);
				fun_s +=  weight->data.vctr[i] * ( log( 
														exp(-(bb->data.vctr[i]))+
														exp((aa->data.vctr[i])-(bb->data.vctr[i])))
													+ (bb->data.vctr[i])) ;
				prob->data.vctr[i] = 1. / (1. + exp(aa->data.vctr[i]) );
				b->data.vctr[i] = -(weighty->data.vctr[i]) *(1 - prob->data.vctr[i]);

				gc += b->data.vctr[i];
			}

			NORMALIZEDMTV(A,b,g,opts,b1);


			slep_dvcopy(x,xp);
			slep_dvcopy(Ax,Axp);
			cp = (*c);
			slep_dvcopy(t,tp);

			for(;;){
				(*c) = sc - gc/L;
				slep_dvcopy(s,u);
				slep_daxpy(-1.0/L, g, u);
				DMEMSET(v->data.vctr,-lambda/L,v->dim,i);
				slep_daxpy(1., s_t, v);

				eppVectorRI(x->data.vctr, t->data.vctr, u->data.vctr, v->data.vctr, ind, x->dim, k);

				slep_dvcopy(x,v);
				slep_daxpy(-1.0,s,v);

				slep_dvcopy(t,v_t);
				slep_daxpy(-1.0, s_t, v_t);
				
				NORMALIZEDMV(A,x,Ax,opts,x1,i);

				slep_dvcopy(Ax, aa);
				slep_xdiagy(y,aa);
				slep_daxpby(-(*c),y,-1.0,aa);

				fun_x = 0.;
				for(i=0; i< aa->dim; i++){
					bb->data.vctr[i] = max(aa->data.vctr[i], 0.);
					fun_x +=  weight->data.vctr[i] * ( log( 
															exp(-(bb->data.vctr[i]))+
															exp((aa->data.vctr[i])-(bb->data.vctr[i])))
															+ (bb->data.vctr[i])) ;
				}


				r_sum = (slep_ddot(v,v) + slep_ddot(v_t,v_t) + ((*c)-sc)*((*c)-sc))/2.0;
				l_sum = fun_x - fun_s - slep_ddot(v,g) - ((*c)-sc)*gc;

				if( r_sum <= 1e-20){
					bFlag=1;
					break;
				}

				if( l_sum <= r_sum * L)
					break;
				else
					L = max(2.0*L, l_sum / r_sum );
			}

			alphap = alpha;
			alpha = (1.0 + sqrt(4*alpha*alpha + 1.0))/2.0;

			valueL->data.vctr[iterStep] = L;

			slep_dvcopy(x,xxp);
			slep_daxpy(-1.0, xp, xxp);
			ccp = (*c) - cp;

			funVal->data.vctr[iterStep] = fun_x + lambda * slep_ddot(t,gWeights);

			if(bFlag)
				break;

			if( opts->tFlag == 0){
				if( iterStep >= 2){
					if( fabs(funVal->data.vctr[iterStep] - funVal->data.vctr[iterStep-1]) < opts->tol)
						break;
				}
			}else if( opts->tFlag == 1){
				if(iterStep>=2){
					if( fabs(funVal->data.vctr[iterStep] - funVal->data.vctr[iterStep-1]) < (opts->tol)*(funVal->data.vctr[iterStep-1]) )
						break;
				}
			}else if(opts->tFlag == 2){
				if( funVal->data.vctr[iterStep] <= opts->tol)
					break;
			}else if(opts->tFlag == 3){
				if(sqrt(slep_ddot(xxp,xxp)+((*c)-cp)*((*c)-cp) + slep_ddot(t,t)+slep_ddot(tp,tp)-2.0*slep_ddot(tp,t) ) <= opts->tol)
					break;
			}else if(opts->tFlag == 4){
				if( sqrt(slep_ddot(xxp,xxp)+((*c)-cp)*((*c)-cp) + slep_ddot(t,t)+slep_ddot(tp,tp)-2.0*slep_ddot(tp,t) )
						<= (opts->tol) * max(sqrt(slep_ddot(xp,xp)+cp*cp + slep_ddot(tp,tp)),1))
					break;
			}else if(opts->tFlag == 5){
				if( iterStep >= opts->maxIter)
					break;
			}
		}
		slepFreeVector(weighty);
		slepFreeVector(xp);
		slepFreeVector(Axp);
		slepFreeVector(xxp);
		slepFreeVector(t);
		slepFreeVector(tp);
		slepFreeVector(v_t);
		slepFreeVector(s);
		slepFreeVector(As);
		slepFreeVector(s_t);
		slepFreeVector(aa);
		slepFreeVector(bb);
		slepFreeVector(g);
		slepFreeVector(prob);
		slepFreeVector(u);
		slepFreeVector(v);

		funVal->dim = iterStep;
		valueL->dim = iterStep;
	}

	/* Reformulated problem + adaptive line search scheme */
	/* Original problem:           min  f(x,c) = - weight_i * log (p_i) + 1/2 rsL2 * ||x||_2^2
	 *                                                      + z * \|x\|_1
	 *
	 * Reformulated problem:   min  f(x,c) = - weight_i * log (p_i) + 1/2 rsL2 * ||x||_2^2
	 *                                                      + z * 1^T t
	 *                                   s.t., |x| <=t
	 **/
	if( opts->mFlag==1 && opts->lFlag==1  && (opts->q > 1.99 && opts->q < 2.01) ){
		int bFlag = 0, iterStep;
		double L = 1.0 /(A->m);
		slep_vector *weighty, *xp, *Axp, *xxp, *aa, *bb, *prob;
		double cp, ccp, gamma, alpha, alphap, beta, gc, sc, cnew;
		slep_vector *t, *tp, *s, *s_t, *As;
		double fun_s, fun_x, l_sum, r_sum, tao;
		slep_vector *u, *v, *g, *xnew, *tnew, *v_t, *Axnew;

		weighty = slepCopyToNewVector(weight);
		slep_xdiagy(y, weighty);

		gamma=1.0;
		
		xp = slepCopyToNewVector(x);
		Axp = slepCopyToNewVector(Ax);
		xxp = slepAllocDenseVector(x->dim);
		slep_zerov(xxp);

		cp = (*c);
		ccp = 0.0;

		t = slepAllocDenseVector(k);
		slep_zerov(t);
		for(i=0; i < k; i++){
			temp = (*x) ; 
			temp.data.vctr = x->data.vctr +(ind[i]);
			temp.dim = ind[i+1] - ind[i];
			t->data.vctr[i] = slep_dnrm2(&temp);
		}
		tp = slepCopyToNewVector(t);

		s = slepAllocDenseVector(x->dim);
		s_t = slepAllocDenseVector(t->dim);
		tnew = slepAllocDenseVector(s_t->dim);
		As = slepAllocDenseVector(Ax->dim);
		u = slepAllocDenseVector(x->dim);
		v = slepAllocDenseVector(x->dim);
		v_t = slepAllocDenseVector(s_t->dim);
		xnew = slepAllocDenseVector(x->dim);
		Axnew = slepAllocDenseVector(Ax->dim);
		g = slepAllocDenseVector(x->dim);
		aa = slepAllocDenseVector(weight->dim);
		bb = slepAllocDenseVector(weight->dim);
		prob = slepAllocDenseVector(weight->dim);

		for( iterStep=0; iterStep< opts->maxIter; iterStep++){
			for(;;){
				if( iterStep !=0){
					alpha = (-gamma + sqrt(gamma*gamma + 4.0*L*gamma)) / (2.0 *L);
					beta = (gamma - gamma * alphap) / (alphap * gamma + alphap * L * alpha);

					slep_dvcopy(x,s); 
					slep_daxpy(beta,xxp,s);

					slep_dvcopy(t,s_t);
					slep_daxpby(-beta, tp, 1.0+beta, s_t);
					sc = (*c) + beta * ccp;

					slep_dvcopy(Ax, As);
					slep_daxpby(-beta, Axp, 1.0+beta, As);
				}else{
					alphap = alpha = (-1.0 + sqrt(5.0)) / 2.0; beta = 0.0;
					slep_dvcopy(x,s);
					slep_dvcopy(t,s_t);
					sc = (*c);
					slep_dvcopy(Ax,As);
				}

				slep_dvcopy(As, aa);
				slep_xdiagy(y,aa);
				slep_daxpby(-sc,y,-1.0,aa);

				gc = 0.;
				fun_s = 0.;
				for(i=0; i< aa->dim; i++){
					bb->data.vctr[i] = max(aa->data.vctr[i], 0.);
					fun_s +=  weight->data.vctr[i] * ( log( 
															exp(-(bb->data.vctr[i]))+
															exp((aa->data.vctr[i])-(bb->data.vctr[i])))
														+ (bb->data.vctr[i])) ;
					prob->data.vctr[i] = 1. / (1. + exp(aa->data.vctr[i]) );
					b->data.vctr[i] = -(weighty->data.vctr[i]) *(1 - prob->data.vctr[i]);

					gc += b->data.vctr[i];
				}

				NORMALIZEDMTV(A,b,g,opts,b1);

				cnew = sc - gc/L;
				slep_dvcopy(s, u);
				slep_daxpy(-1.0/L, g, u);
				DMEMSET(v->data.vctr,-lambda/L,v->dim,i);
				slep_daxpy(1.0, s_t, v);

				eppVectorRI(xnew->data.vctr, tnew->data.vctr, u->data.vctr, v->data.vctr, ind, xnew->dim, k);

				slep_dvcopy(xnew,v);
				slep_daxpy(-1.0,s,v);

				slep_dvcopy(tnew,v_t);
				slep_daxpy(-1.0,s_t,v_t);

				NORMALIZEDMV(A,xnew,Axnew,opts,x1,i);

				slep_dvcopy(Axnew, aa);
				slep_xdiagy(y,aa);
				slep_daxpby(-cnew,y,-1.0,aa);

				fun_x = 0.;
				for(i=0; i< aa->dim; i++){
					bb->data.vctr[i] = max(aa->data.vctr[i], 0.);
					fun_x +=  weight->data.vctr[i] * ( log( 
															exp(-(bb->data.vctr[i]))+
															exp((aa->data.vctr[i])-(bb->data.vctr[i])))
															+ (bb->data.vctr[i])) ;
				}

				r_sum = (slep_ddot(v,v) + slep_ddot(v_t,v_t) +(cnew-sc)*(cnew-sc))/2.0;
				l_sum = fun_x - fun_s - slep_ddot(v,g) - (cnew-sc) * gc;

				if(r_sum <= 1e-20){
					bFlag=1;
					break;
				}

				if(l_sum <= r_sum * L)
					break;
				else
					L= max(2.0*L, l_sum / r_sum );

			}// End for(;;)

			gamma = L * alpha * alpha;
			alphap = alpha;

			valueL->data.vctr[iterStep] = L;

			tao = L * r_sum / l_sum ;
			if(tao >= 5.0)
				L*= 0.8;

			slep_dvcopy(x,xp);
			slep_dvcopy(xnew, x);
			slep_dvcopy(x,xxp);
			slep_daxpy(-1.0, xp, xxp);
			slep_dvcopy(Ax, Axp);
			slep_dvcopy(Axnew, Ax);
			slep_dvcopy(t,tp);
			slep_dvcopy(tnew, t);
			cp = (*c);
			(*c) = cnew;
			ccp = (*c)-cp;

			funVal->data.vctr[iterStep] = fun_x + lambda * slep_ddot(t,gWeights);

			if(bFlag)
				break;

			if( opts->tFlag == 0){
				if( iterStep >= 2){
					if( fabs(funVal->data.vctr[iterStep] - funVal->data.vctr[iterStep-1]) < opts->tol)
						break;
				}
			}else if( opts->tFlag == 1){
				if(iterStep>=2){
					if( fabs(funVal->data.vctr[iterStep] - funVal->data.vctr[iterStep-1]) < (opts->tol)*(funVal->data.vctr[iterStep-1]) )
						break;
				}
			}else if(opts->tFlag == 2){
				if( funVal->data.vctr[iterStep] <= opts->tol)
					break;
			}else if(opts->tFlag == 3){
				if(sqrt(slep_ddot(xxp,xxp)+((*c)-cp)*((*c)-cp) + slep_ddot(t,t)+slep_ddot(tp,tp)-2.0*slep_ddot(tp,t) ) <= opts->tol)
					break;
			}else if(opts->tFlag == 4){
				if( sqrt(slep_ddot(xxp,xxp)+((*c)-cp)*((*c)-cp) + slep_ddot(t,t)+slep_ddot(tp,tp)-2.0*slep_ddot(tp,t) )
						<= (opts->tol) * max(sqrt(slep_ddot(xp,xp)+cp*cp + slep_ddot(tp,tp)),1))
					break;
			}else if(opts->tFlag == 5){
				if( iterStep >= opts->maxIter)
					break;
			}
		}// End for(iterStep)

		slepFreeVector(weighty);
		slepFreeVector(xp);
		slepFreeVector(Axp);
		slepFreeVector(xxp);
		slepFreeVector(aa);
		slepFreeVector(bb);
		slepFreeVector(prob);
		slepFreeVector(t);
		slepFreeVector(tp);
		slepFreeVector(s);
		slepFreeVector(s_t);
		slepFreeVector(As);
		slepFreeVector(u);
		slepFreeVector(v);
		slepFreeVector(g);
		slepFreeVector(xnew);
		slepFreeVector(tnew);
		slepFreeVector(v_t);
		slepFreeVector(Axnew);

		funVal->dim = iterStep;
		valueL->dim = iterStep;
	}// End if

	slepFreeVector(weight);
	slepFreeVector(b);
	slepFreeVector(x1);
	slepFreeVector(Ax);
	slepFreeVector(ATb);
	slepFreeVector(b1);
	if( ! (hasGot(opts, SlepOPT_gWeights)) )
		slepFreeVector(gWeights);
	return;
}